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ABSTRACT 



In a recent paper (Pearson, Horne & Skidmore 
2003), we explained the unusual flaring activity 
of the cataclysmic variable (CV) system AE Aqr 
in terms of the aftermath of the collision between 
two gas clouds. We modelled the resulting fireball 
numerically, comparing the results to analytic ap- 
proximations for the optical lightcurves and con- 
tinuum spectra and to observed lightcurves and 
spectra. We set out here an improved analytic 
calculation for an expanding fireball with an LTE 
ionization structure and compare our results to 
more detailed numerical simulations. We fit our 
results to multiwavelength lightcurves of two very 
different systems; deriving values for the physical 
conditions in them. 



creting systems and, in particular, those where an 
accretion disk is present. This flickering process 
may well be associated with the anomalously high 
viscosity present in these systems and represent 
the effect of magnetic reconnections such as that 
from the viscosity mechanisms of Hawley & Bal- 
bus (1991). Such a sudden localised deposition 
of energy would raise the gas temperature and the 
consequent overpressure would give rise to a lo- 
cal expansion of the disk material. With material 
cooling adiabatically, this expansion would rapidly 
become supersonic. 



Flickering and flaring occurs across the whole 
range of astrophysical systems from stars to active 



Observations by Bruch (2000) suggest that 
the flickering in CVs is not uniformly distributed 
across the disk but instead is associated with 
the stream-disk impact point and the innermost 
boundary layer region of the disk. Similarly Pat- 
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terson (1981) found that the flickering in HT Cas 
was associated with the inner part of the accretion 
disk. These findings, however, are contradicted by 
Welsh & Wood (1995) for HT Cas, by Baptista 
(2004) for V2051 Oph and Baptista (2002) for the 
Low-Mass X-ray Binary (LMXB) X1822-371. All 
these studies show flickering to arise from a range 
of disk radii. 

The terms "flickering" and "flaring" are often 
used interc;hangably when describing the stochas- 
tic variability of accreting sources and, even when 
defined in a particular field, are not necessarily 
used consistently. Warner (1995) and Baptista 
(2004) describe flickering in CVs as the continu- 
ous, random brightness fluctuations of 0.01-1 mag 
on timcscales of seconds to dozens of minutes al- 
though the exact numerical ranges differ between 
the two. In contrast the unusal CV AE Aqr is 
universally described as a "flaring" source since it 
has rarer by brighter events that often occur in 
batches. These events can raise the total optical 
luminosity of the system by factors of 2-3 contrast- 
ing with the 5-20% typical of CV flickering (Bruch 
1992) but on a similar timescale. Consequently, we 
adopt the convention of describing the small am- 
plitude, continuous variations as "flickering" and 
reserve the term "flaring" for larger scale events. 

The underlying model used to reproduce the 
AE Aqr lightcurves was based on consideration 
of the flux emitted by a hot, spherically symmet- 
ric ball of gas as it expanded and cooled. While 
we saw that the "expanding fireball" situation in 
AE Aqr could arise from coUisional heating of gas 
blobs, the same situation might arise from the 
very different causes oiitlined above. As a result, 
these "fireball" models may have a much wider 
range of applicablility than simply the unusual be- 
haviour of AE Aqr. Accordingly, in this paper, we 
will compare analytic expressions derived for the 
lightcurves of such expansions to observations of 
SS Cygni and SN 1987A. 

SS Cyg is a member of the dwarf nova subclass 
of cataclysmic variables. It consists of a 1.2Mq 
white dwarf accreting material through an accre- 
tion disk from a OJIMq K5V secondary that loses 
material via Roche lobe overflow (Ritter & Kolb 
2003). SS Cyg was the second member of this 
subclass identified (Wells 1896) and has had a 
near continuously monitored lightcurve for over a 
century (Warner 1995). Like all dwarf novae, 



SS Cyg exhibits optical variability over a large 
range of timescales e.g. outbursts (Am « 3.5 mag, 
f ~ 40 days), orbital modulations (Am w 0.5 mag, 
P = 0.275130 days), flickering (Am 0.01 - 0.2 
mag, t ~ 1 min) and Dwarf Nova Oscillations 
(Am « 0.02 mag, t ^ 2.8 - 10.9 s) (Mauche & 
Robinson 2001; van Teeseling 1997; Honey et al. 
1989; Warner 1995, 2004). We shaU concentrate 
our modelling on the lightcurve from a flickering 
event. 

SN 1987A was the brightest supernova observed 
since Kepler's in 1604. The lightcurve can be bro- 
ken into three phases: a "flash" lasting a few 
hours, a subsequent "bump" lasting around 4 
months and a final exponential "tail" powered by 
radioactive decay of unstable nuclei in the ejecta. 
It was a type II supernova with an identified pro- 
genitor (Sk-69''202) that was a blue supergiant. 
It is believed that it lasted ~ lO'' y as a 20-30Mq 
main sequence blue supergiant star, before becom- 
ing a red supergiant for ~ 10^ y during which 
phase it lost 3-6M0 before finally becoming a blue 
supergiant again for a few thousand years (Mc- 
Cray 1993). The rebrightening "bump" phases 
arises from a photosphere moving out with ex- 
pelled material and eventually reversing as the 
material becomes optically thin again. It is in- 
teresting to note that the photosphere maintained 
a roughly constant temperature ~ 6000 K close 
to the recombination temperature for hydrogen 
at the expected densities in the expanding shell. 
The models of AE Aqr flares in Pearson, Horne 
& Skidmore (2003) showed a similar isothermal 
behaviour. 

2. Model Assumptions 

We briefly recap here the dynamical model de- 
veloped in Pearson, Horne & Skidmore (2003) but 
the interested reader is referred there for a more 
detailed derivation. Let us assume a spherically 
symmetric expansion of a Gaussian density profile 
with radial velocity proportional to the distance 
from the centre of the expansion. We define 



as a dimensionless measure of the radius r in terms 
of the current lengthscale of the Gaussian a. The 
expansion factor (3 acts as a dimensionless time, 
being the constant of proportionality between the 
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current and fiducial scalelength ao- Hence, 

(2) 

ao 

consistent with our above assumptions regarding 
velocity and implying 



f3 = l + H {t-to) 



(3) 



where to is the time at which all the fiducial values 
are determined and H is an "expansion constant" 
setting the speed of the expansion. For simplic- 
ity, we consider only the case of uniform 3-D ex- 
pansion, avoiding factors such as the angle of the 
observer to lines of symmetry. We also restrict 
ourselves to a spatially uniform temperature dis- 
tribution since, for example, a power law distribu- 
tion gives a result for the later integration in (22) 
in terms of the hypergeometric function 2-Fi- We 
assume a power law with index F for the temporal 
dependance of temperature. F = corresponds to 
an isothermal expansion and F = 2 to the adia- 
batic case. 

In summary then, we have 

-r 



T = Top- , 
P = Po li'^ 



where 
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(4) 
(5) 
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and M is the total mass involved in the expansion. 

Using v{r) = Hro = Hrj3~^ we can integrate 
\pv'^ over all space to derive the total kinetic en- 
ergy of the expansion 
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defining vq = Huq, the speed of expansion at r = 
a. 

3. Theoretical Lightcurves and Spectral 
Distributions 

The radiative transfer equation has a formal so- 
lution under conditions of LTE 
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dr, 
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where / is the intensity of the emerging radiation, 
B is the Planck function and r is the optical depth 
measured along the line of sight from the observer. 
For cases where the source function is everywhere 
the same this becomes 



I = B (1-e-^). 



(11) 



We define x as the distance from the fireball 
centre toward the observer, and y as the distance 
perpendicular to the line of sight. The above line 
integral (equation (11)) gives the intensity I{y) for 
lines of sight with different impact parameters y. 
The fireball flux, obtained by summing intensities 
weighted by the solid angles of annuli on the sky, 
is then 

fiX)=J^^I{y)^dy (12) 

where d is the source distance. 

We can calculate the evolution of the continuum 
lightcurve using expressions for the linear absorp- 
tion coefficient. The free-free absorption coefB- 
cient (per unit distance) can be written as 



Kff = Kq 
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where 



Ko = 3.692 X 10-2 Q2 (sj) 



(13) 
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(Keady & Kilcrease 2000), gs is the free-free 
Gaunt factor (Gaunt 1930) and Q the ionic 
charge. We wish to retain the explicit frequency 
dependence but otherwise follow a parallel deriva- 
tion as for Kramers' opacity. For a mixed ele- 
mental composition, then, we sum over all species 
(assumed fully ionized) 



all ions 



(1 - Z) -gs 



(15) 



(Bowers & Deeming 1984). Here ga is a mean 
Gaunt factor (close to unity) and Z the metal mass 
fraction. We also have 



P 



p{i + x) 

2 mn 
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(Bowers & Deeming 1984) where X is the hy- 
drogen mass fraction. Combining and writing the 
correction for stimulated emission as 



= l-e-(fe) 



(17) 
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(18) 



gives 

where we have defined 

Ki = 6.695 X 10-^1 (1 - Z) (1 + X)gs (SI). (19) 

Inserting our density profile (equation (5)) we ar- 
rive at the result 

Kff = — r ^ 6-2" . 20 



On dimensional grounds if on no other, a sim- 
ilar expression to equation (18) and hence (20) 
must also exist for bound-free opacity when it 
dominates (when most species are fully recom- 
bined). Textbook derivations of Kramer's opacity 
(eg. Bowers & Deeming 1984) show us 



«;i,bf (xZ {1 + X) gu 



(21) 



where gbf is the bound-free Gaunt factor (Gaunt 
1930). Useful tables for gbf have been calculated 
by Glasco & Zirin (1964). The constant of propor- 
tionality, however, is more difficult to determine 
than for free-free, since it must account for the 
ionization edges in the absorption and, in partic- 
ular, the change of energy level populations with 
temperature for each ion. We shall assume a rela- 
tion analagous to equation (18) also holds for the 
situation of mixed ionized and recombined species 
where both forms of opacity contribute. 

The optical depth parallel to the observers line 
of sight 
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where the optical depth on the line of sight 
through the centre of the fireball is 



TO 
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(25) 



and the time at which the fireball becomes opti- 
cally thin along this line of sight is 



/3c 



(26) 



It should be noted that in pulling Ki(/ii, ji^) out of 
the integral (23) we have implicitly assumed that 
the spatial variation of the ionization fraction has 
negligible impact on the behaviour. This is an 
assumption that we shall return to later. 

In Pearson, Home & Skidmore (2003) we ap- 
proximated the fiux received from the fireball by 
splitting it into two components as viewed on the 
plane of the sky: an optically thick central region 
bounded hy y = y^a (see equation 44) and an op- 
tically thin surrounding. Here, however, we cal- 
culate the integral exactly. Combining equations 
(11), (12) and (24) with a change of integration 
variable to u = T(y), we have 
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Expressing the exponential in its series form we 
integrate to find 
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ra^B^ (-l)("+i)ro" 

1^ ^ (28) 



which can also be rewritten in the form of standard 
functions 



/ 



2 d? 



-(£;i(ro)+7 + ln(ro)) (29) 



where Ei is the first order exponential integral 
and 7 « 0.577216 is Euler's gamma constant 
(Abramowitz & Stegun 1972; Jefi'reys & Jefi'reys 
1956; Press et al. 1992). This has the form of 



/ = S 5(to), 



(30) 



where Q. = is the solid angle subtended by the 
current standard deviation (a/v^) of the Gaussian 
density profile and the "Saturation Function" 



5(ro) = (i;i(ro)+7 + ln(To)) 



(31) 



is plotted in Figure 1. We note the asymptotic 
limits 
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for r «; 1 
for T > 1. 



(32) 



The intensity of emitted radiation as a function of 
impact parameter is plotted for different times in 
Figure 2. We can see how the flux saturates at 
the black body function for highly optically thick 
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impact parameters. In comoving coordinates this 
region gradually shrinks and at later times the en- 
tire emission becomes optically thin. 

The total flux given by equation (29) is plotted 
against time in Figure 3 for different values of F, 
using Pc = 1. Figure 4 shows lightcurves at differ- 
ent wavelengths assuming F = (isothermal) and 

= (A/5000 A)i. The effect of the expanding 
photosphere can be seen from the optically thick 
contribution to the total flux, the latter of which 
rises to a peak at /3pk. Initially the optic;ally thick 
flux is the dominant source and the emission rises 
as the emitting area grows while the photosphere 
is advected outwards. Eventually, tlw^ decreasing 
density causes the photosphere to reverse and be- 
gin to shrink, reaching zero size at /3 = /3c- The re- 
maining emission comes from the optically thin re- 
gion surrounding the optically thick circle as seen 
on the plane of the sky. This optically thin emis- 
sion dominates at late times. 

Recalling that a = /3ao, for the special case of 
F = (isothermal expansion) we can neglect the 
time variation of the Planck function and differen- 
tiate (29) with respect to /3 to find a condition for 
the maximum flux 

4(i;i(To)+7+ln(To))-(10-r)(l-e^^«) =0. (33) 

This can be solved numerically to find To,pk = 
6.8204. From (25) we then find /3pk = 0.6811/3c. 
We can also differentiate (29) with respect to /? for 
F ^ although the form is much more untidy and 
we must solve for /3 directly. We plot /3pk against 
F in Figure 5 continuing to neglect any time de- 
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Fig. 1. — The Saturation Function S. 
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Fig. 2. — Comparison of the intensity profiles at 
different times calculated using (24) and (11) for 
a F = (isothermal) evolution. 
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Fig. 3. — The temporal behaviour of the to- 
tal fiux at 5000 A for different cooling indices 
F = —2,-1,0,1,2, assuming = 1 and Tq = 
15 000 K. The y-axis is parameterised in multiples 
of /o = If^ 5^(5000 A, 15 000 K). 
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pcndcncc of e. Wc plot /3pk against wavelength in 
Figure 6 for several values of F. 

Bruch (1992) found a correlation between the 
amplitude of flickering in different wavebands for 
several CVs (see inter alia his Figure 2). Such 
a correlation could arise from an isothermal evo- 
lution at a consistent temperature. This can be 
understood from (29). Since the peak fltix occurs 
at the same central optical depth independent of 
wavelength, then for two wavelengths Ai and A2 
the peak fluxes are related by 

/pk,l ^ / /^pk.l V Bx,l ^g^^ 
/pk,2 V/^pk,2/ 

(,^0,2) (A2) (e^-i) ^^^^ 

2 

(36) 

where C2 = hc/k and assuming no complications 
such cLS cL Balmcr edge between them. Wc plot 
the predicted values for different temperatures as 
stars: alongside the data taken from Tabic 4 of 
Bruch (1992), in a color-color diagram in Fig- 
ure 7 (nb. we used B\ rather than in the 
above derivation purely for comparison with this 
dataset). 

This figure shows that a simple application of 
our model enables us to reproduce the observed 
range of V-B reasonably well using T as a free 
paramemter. The U-B colour however, is gener- 
ally underprcdicted. This results from not mak- 
ing allowance for the increased opacity above the 
Balmer jump when evaluating /3c in deriving 36. 
The size of the Balmcr jump depends upon a 
non-trivial combination of mass, lengthscale and 
temperature to produce an effective value of ki. 
Any individual system may produce flickers with 
a range of physical parameters causing it to pro- 
duce points scattered across this diagram. 

Using numerical methods outlined later, we de- 
rived values for ki for several temperatures over a 
range of densities. These are plotted in Figure 8. 
We can see how, for a large range of temperatures 
of interest and for densities ranging over several or- 
ders of magnitude, ki is well represented by a con- 
stant value. At lower temperatures the free-free 
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Fig. 4. — The temporal behaviour of the total flux 

3 

at several wavelengths assuming f3c = ( — 1 , 

V5000 A/ 

To = 15 000 K and F = 0. The y-axis is again 
parameterised in multiples of /q. The optically 
thick contribution at 5000 A is plotted as a dashed 
line and the optically thin contribution as a dot- 
dashed line. 
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Fig. 5. — The time of peak flux at 5000 A plotted 
against F, assuming (3c = 1 and Tq = 15 000 K. 
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Fig. 6. — The time of peak flux plotted against 

wavelength for various values of F assuming To = 
15 000 K. The values are normalised to the peak 
time at 6000 A in each case. 
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Fig. 7. — The amplitude ratio (a U-B v V-B color- 
color diagram) for flickering in several CVs ob- 
served by Bruch (1992). The predicted values arc 
plotted as stars for different temperatures ranging 
from 8000 K (extreme right) to 24 000 K in 2000 K 
steps. 



opacity transitions to bound-free opacity at lower 
densities than it does at higher temperatures. At 
these lower temperatures we could introduce sig- 
nificant errors if we estimate ki from a position 
that has one form of opacity when the other is in 
fact dominant. It is very unlikely that this would 
be the case, however, if wc estimate the ionization 
state from conditions at a suitable position and 
consider how Ki smoothly changes with the ion- 
ization. Qualitatively then, at the very least, wc 
can be confident of the ability of the analysis to 
allow us to examine the behaviour of the proposed 
mechanism and to predict fluxes to within factors 
of order unity. 

4. Analysis 

Equation (29) gives us the complete temporal 
and frequency behaviour for the emergent flux in 
terms of tq. We must use the Planck function 
appropriate to the time-dependent behaviour of 
the temperature from (4). Equation (25) gives us 
To as a function of /3 and /3c. The behaviour of the 
flux thus rests on the character of the parameter 
/3c- We know functionally that 

/3c = /3c(Ki,e,!/-3), (37) 
Ki = Ki{r], 1^,(3), (38) 
e = e{,^,(3) (39) 

and we shall examine each of these dependencies 
and appropriate levels of approximation in each 
case. 

4.1. Global Behaviour /3c = const 

In examining the behaviour of the flux as a 
function of time it is instructive to restrict our- 
selves initially to a single wavelength and ignore 
spatial or time dependency in /3c. 

We see from equations (25) and (28) that the 
late-time behaviour {kT <C hu) of the lightcurve is 
determined by the temperature index F such that 
/ (3~^ or (3~^ for isothermal or adiabatic cooling 
respectively. 

4.2. Spatial Dependence (3c = (3c{ri) 

4-2.1. One Zone Model 

We know depends on ki which in turn de- 
pends on the ionization structure across the ex- 
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pansion. At the outset we acknowledge that a de- 
tailed solution accounting for the changing ioniza- 
tion fraction across the density profile would re- 
quire numerical integration such as that reported 
in Pearson, Horne & Skidmore (2003) for AE Aqr. 
We are looking for simpler, more easily calculable 
approximations that avoid such detailed methods. 
To this end we do not include here any explicit al- 
lowance for the variation of ionization states across 
the profile; instead we approximate the behaviour 
by calculating the conditions at a suitable point 
within the expansion. With a Gaussian density 
profile we might expect the flux to be dominated 
by the region close to the photosphere. As a result 
we could approximate the ionization structure by 
the solution to the Saha equation under the con- 
ditions prevailing there. 

To find the position of the observed photo- 
sphere we need to integrate along a given line of 
sight until optical depth unity. From our earlier 
derivation we can see the optical depth down to 
any height x' is 

r(,)=ro^e-(*)YJe-(t)^dQ (40) 

which, setting r(y) = 1 for the special 
gives 

e^il? = erfc (v^^) • (41) 

This defines the locus of points on the photo- 
spheric surface. Since the largest contribution to 
the opacity integral comes from the region of high- 
est density, we will introduce least error by eval- 
uating Ki there. This will clearly occur along the 
central line of sight y = and hence we arrive at 



erfc(A/2r?ph) = — 
To 



(42) 



which we must solve numerically. However, we 
must remember that it is possible for the pho- 
tospheric surface to lie behind the density peak. 
This will just occur when t(0) = 1 at 77 = 0, im- 
plying /3 = 2 (10^) /3c. Ultimately then, we have 



Peval 



po e-''p'> for /? < 2^-^)^^ 
po otherwise. 

(43) 

It must be emphasized that we are using this 
position solely to evaluate a typical ionization 



state and, hence, to find a suitable approximate 
value for k\. The effect of the density profile on 
the optical depth is accounted for explicitly in the 
derivations in section 3. 

However, an alternative 'typical' place to es- 
timate K\ is at the point of maximum density 
(x = 0) along the limiting impact parameter y-^ 
that remains optically thick. We can derive t/m by 
setting (24) equal to unity, giving 



2/m = do /? 



In To 



(44) 



where Tq > 1 ie. for times P < Pc- Hence, 



(-2-r) 



Peval 



1 Po 



(lo-r) 



for f3<f3c 
otherwise. 



(45) 



4-2.2. Pure hydrogen case - semi- analytic solu- 
tion 

For simplicity, let us consider a fireball with 

almost pure hydrogen composition (nb. if Z = 
then equation (21) implies Ati^bf = 0), with an 
ionization fraction 



(46) 



Ignoring any contribution to opacity from H~ and 
noting K (X p^, the total opacity coefficient be- 
comes 

Kl « Kl.flF + {1- if Kl,bf • (47) 

For hydrogen, the Saha equation can be approx- 
imated by 



Uj Tie 

rin 



n — n-i 



2.4 X 10^^ Ti exp(-1.58 x 10^ T'^) m" 



A{T) 



(48) 
(49) 



(Lang 1980) where we have tidied the RHS into 
the function A{T). Hence, we derive the quadratic 



nf + A{T) rii - A{T) n = 



(50) 



and by restricting ourselves to the positive root 
arrive at 



(51) 
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Wc sec then that the ionization fraction l and 
parameter /3c are interdependant and should be 
solved for iteratively to ensure consistency. Thus 
wc can determine the lightcurve behaviour by ini- 
tially calculating a value for /3c from (26) with a 
judiciously chosen value for k\ . Thereafter we use 
the number density from (43) to calculate the ion- 
ization fraction from (51). This then allows us to 
find Ki from (47) and hence /3c again from (26) to 
better accuracy. We can repeat this iteratively to 
achieve the desired accuracy. This value can then 
be used to calculate the lightcurve behaviour at a 
given wavelelength. 

4.-2.3. Mixed composition 

For the more realistic case of mixed composi- 
tion, wc cannot calculate the ionization state at a 
given position analytically. Instead we must iter- 
atively solve the network of Saha equations under 
the prevailing physical conditions to enable us to 
determine ki . This solution can then iterated with 
/3c around the equations (26), (43) and ki loop in 
a similar way to the above. 

4.2.4. Three Zone Model 

The above method works well in situations 

where either free-free or bound-free opacity is 
dominating. Unfortunately, as we noted at the 
end of section 3, in the situation of mixed opacity 
sources we can introduce significant errors if we 
use a value for ki assuming the incorrect source. 
We can rescue much of the above formalism, how- 
ever, if we split the spatial profile into 3 zones. 
We are almost bound to improve our calculation 
regardless of where we place the boundaries of 
the zones but clearly it makes sense to try to en- 
sure that we have free-free dominated (outer re- 
gions), bound- free dominated (inner regions) and 
mixed (intermediate regions) zones. Assuming 
that hydrogen species provide the dominant opac- 
ity source, wc select the boundary between these 
regions by the hydrogen ionization fractions L = 
0.1 and t = 0.9 which occur at 770.1 and 770.9. 

Specifically then, we need to replace the inte- 
gral in equation (23) with one over the three types 
of zone and rework equation (26) to redefine /3c. 
Calculating a;o.i and 0:0.9 using y = or j/m as 



desired, the integral in (23) thus becomes 



poo 








J —oc 
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f-OC 




'xo.9 



= 2/^i,ff erfc {V2^) 

+2Kl,m 
+2Kl,bf 



erfc (72^) -erfc (72^) 
"l_erfc(V2^) 



which we use to replace a factor m \^Tr/2 in equa- 
tion (26). 

The question remains of how to rapidly find 
values for 770.1 and 779.9. Again, in the spirit of 
finding an easily calculable approximation, and 
noting that even inaccurately determining these 
boundaries will still improve our integral calcula- 
tion, let us consider a situation where all species 
other than hydrogen remain fully ionized. Using 
standard methods, we derive 



Ate 



1 -I- X (2(. - 1) ' 



(54) 



Incorporating this into equation (48) along with 
our density profile, we can show 



2 777,^ 



[l-|-(2 i-1) X] 



-A{T) (55) 



which we can rearrange and solve directly for 77^. 
We note the particularly simple form this reduces 
to for L = 0.5. 

For similar reasons to the single zone model, 
we evaluate Kbf and Kg at the point of highest 
density in their region (eg. either at 770.9 or 77max 
for free- free) . Given the rapid density variation in 
the mixed zone, we evaluate at 770.5 or 7/max as 
appropriate. The possible integration schemes are 
illustrated schematically in Figure 9. 

4.3. Frequency Dependence /3c = /3c(?7, J^) 

The parameter /3c has a direct frequency de- 
pendence from the in opacity and also depen- 
dence through both ki and e. For regions where 
the Gaunt factors are slowly changing functions 



(52) 



(53) 
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of frequency eg. along the Paschen continuum in 
the 4000-8000 A range or when free-free opacity is 
dominating, we ignore the small error introduced 
by neglecting the contribution of gff and gbf to ki. 
Instead we need only correct for the direct and e 
frequency dependence of /3c- Thus, 

/3cM«/3cM(^)"^(^)"^. (56) 

With this correction we have only to calculate /3c 
accurately at a single wavelength with the iterative 
method outlined above. The lightcurve behaviour 
as a function of both wavelength and time then 
follows immediately from the combination of (56) 
and (29). 

4.4. Time Dependence /3c = /3c(?7, z^, /3) 

The dependence of /3c on time again comes 
through both ki and e (if F 7^ 0). In principle, we 
can use a similar expression to equation (56) also 
to correct for the e time-dependence. However, 
the exponential nature of the expression renders 
the form of the lightcurves more complex and less 
instructive than that derived above. As a result, 
this correction is probably best included only in 
numerical solutions. More seriously, we cannot, in 
general, predict the future ionization structure for 
a mixed composition gas from its current state. 
Thus, the time-dependence of ki can, in general, 
only be included with numerical solution at a se- 
ries of different times /3. The exception to this rule 
is when we can be sure that the dominant opacity 
is and will remain free-free throughout the time 
considered. In this case ki is a constant in time 
and we can use the same rapid approximation as 
the previous section. 

4.5. Comparison of Integral Methods 

We can lift the restriction to the purely free- free 
opacity case by iterative numerical solution of the 
Saha equations for a gas of mixed composition at 
each time in the lightcurve. Prom the ionization 
profile of each species we can in principle deter- 
mine the opacity at any point. 

We compared the results achieved by the two 
integration lines of sight y ~ Um or and using 
1- or 3-zone integration schemes. The total con- 
tinuum opacity was calculated numerically using 



using the methods of Gray (1976) with the excep- 
tion of H~ bound-free (Geltman 1962) and free- 
free (Stilley & Callaway 1970), He" (McDowell, 
Williamson & Myerscough 1966) and Hel bound- 
free (Huang 1948). Calculating /3c from this opac- 
ity we can iterate to a consistent solution at each 
time. Altering the number of zones used yielded 
virtually no discernable effect on the predicted 
lightcurves. For temperatures around 16 000 K, 
the two lines of sight considered produced results 
differing by at most around 0.5%. However, at 
lower temperatures the predicted fluxes could dif- 
fer more, reaching around 5% at 10 000 K (see 
Figures 10 and 11). 

5. Comparison to Observation - Deducing 
Fireball Parameters 

The model was fitted to lightcurves for the fiick- 
ering in the dwarf nova SS Cyg and the 'bump 
phase' of SN 1987A. We employed a minimiza- 
tion amoeba code (Press et al. 1992) to derive the 
best-fit values for M, uq, Tq, to, H, and T. We use 
the approximation that the ionization fractions for 
all the species are well represented by the condi- 
tions at X = 0,2/ = j/m and use the single zone 
integration approach to arrive at a self-consistent 
solution for /3c and ki there at each time. 

5.1. SS Cygni 

Rapid spectroscopy of SS Cyg was carried out 
using the Low Resolution Imaging Spectrograph 
(LRIS: Oke et al. 1995) on the 10-m Keck II tele- 
scope on Mauna Kea, Hawaii, between UT 09:38 
and 09:56 on 6 July 1998, covering orbital phases 
0.8276 to 0.8715. The instrumental set and data 
reduction were the same as described by Steeglis 
etal. (2001), O'Brien etal. (2001) and Skidmore 
et al. (2003). 14,309 spectra were obtained us- 
ing the rapid data acquisition system. The spectra 
covered 3259-7985 A with 2.4 A pixel"^ dispersion 
and a mean integration time of 72.075 ms and no 
dead-time between individual spectra. Further de- 
tails of these observations are given in Skidmore 
et al. (in prep). A particular flickering event was 
isolated between UT 9:39:10 and UT 9:43:44 and 
flux from other sources removed by fitting, for each 
wavelength, a low order polynomial to the fluxes 
both before and after the event. Lightcurves were 
formed from the mean flux in the regions 3590- 
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3650 A, 4165-4270 A, 4520-4620 A, 5100-5800 A, 
5970-6500 A and 7120-7550 A at 2 s resolution. 

Model lightcurves were computed using a dis- 
tance of 166 pc (Harrison et al. 1999) and assum- 
ing solar composition. It appears that the dis- 
agreement is dominated by the systematic error of 
the oversimplicity of our model rather than obser- 
vational errors. Accordingly, we carried out fits 
both weighting the data points according to their 
formal observational errors and with equal weight- 
ing. We plot the derived analytic lightcurves 
alongside the observations in Figures 12 and 13. 
From each set of derived parameters, we carried 
out a detailed numerical integration of the opacity 
calculated with an LTE ionization varying across 
the density profile as described in Pearson, Horne 
& Skidmore (2003). These models allowed us 
to generate a timeseries of continuum spectra and 
"numerical" lightcurves for comparison to the data 
in each case. These numerical lightcurves are also 
plotted in the Figures. The fitted value of T is 
close to zero in both cases and so wc refitted the 
data fixing P = exactly. These results are shown 
in Figures 14 and 15. The best-fit parameters for 
each case are summarised in Table 5.1. 

Although we would expect the numerical 
lightcurves to more accurately represent reality, 
the analytic lightcurves generally fit the data bet- 
ter. This is unsuprising since the parameters used 
in both cases have been optimized for the ana- 
lytic forms. The numerical lightcurves do seem 
closer to the data in the 3615 Awindow in Fig- 
ures 12 and 13. However, Figures 14 and 15 and 
our other unpublished lightcurves suggest that 
this is serendipitous. The difficulty in reproduc- 
ing the Balmor jump may well reflect that wc have 
5 fitting points on the Paschcn Continuum (wave- 
lengths longer than 3646 A) and only one above 
the Balmer jump. Since the continuum level will 
have a close to relationship between discon- 
tinuities, fixing the level at several points may 
well bias the parameters towards an accurate fit 
here rather than the more complex interplay of 
parameters required to accurately reproduce the 
Balmer jump. In ideal circumstances we would 
like several more points at shorter wavelengths to 
address this issue. 

The derived masses, M ~ 1.6 x 10^^ kg are 
equivalent to 400 s of the estimated mean 
mass transfer rate from the secondary of M w 



4 X 10" kg s-i and - 2 X 10"* of the total 
disk mass Mdisk « 7 x 10^° kg (Schreiber & 
Gansicke 2002). Similarly, the lengthscale of 
the expanding region, a ^ 9.1 x 10^ m is ^1.5% 
of the estimated disk radius i?disk ~ 6 x 10* m 
(Schreiber & Gansicke 2002). The kinetic en- 
ergy is much greater than the thermal energy 
i?th « 4 X 102'^ (js^) (20WK) J but com- 
parable to that of a putative, moderately strong, 
magnetic field in a sphere of radius oq, Ej^ag = 

6 X lO^Mw^)' (e^f J- 

We generated full optical spectra from the '6- 
parameter weighted' set of values and compare 
them to the observed spectra at three differ- 
ent times in Figure 16. The model data have 
been convolved with the instrumental blurring of 
9.8 A FWHM. The results show how parameters 
derived from the analytic forms for the continuum 
lightcurves can subsequently be used to repro- 
duce spectra. The agreement with observation at 
early times is remarkable. At the peak flux and 
later, however, the models predict more metal 
and stronger lines than arc observed, particularly 
in the 3000-5000 A range. Many of the relative 
line strengths, however, are still reproduced. 

5.2. SN 1987A 

Given the similarities in the mechanisms un- 
derlying the flickering lightcurves in our mod- 
els and those of supernovae, we attempted to fit 
the UBVRI photometric data for SN 1987A ob- 
tained from Catchpole et al. (1987). The data 
were converted to fluxes using the standard values 
and effective wavelengths for the bands in Mur- 
din (2001). The fits assumed an interstellar red- 
dening E(B-V)=0.15 (Hamuy et al. 1988; Fitz- 
patrick 1988) and assumed the same distance to 
the LMC of 50.1 kpc as Catchpole et al. (1987). 
Given the results of sophisticated NLTE spectral 
and lightcurve modelling (eg. Mitchell et al. 2002) 
we should not expect an accurate match between 
our simple models and the SN 1987A observations. 
However, it will allow us to highlight the under- 
lying similarities of these problems and to gauge 
the robustness of our derived parameters. The fits 
used only the data after JD 2446875.0 to elim- 
inate the flash phase and the data were all as- 
signed equal weights. For the sake of continuity 
with the SS Cyg fit, we continued to allow 6 fit 
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Model 


M 


ao 


To 


to 


H 


r 


vo 


-Ekin 


SS Cyg 


(10^^ kg) 


(10" m) 


(10* K) 


(UT) 


(s-^) 




(km s-i) 


(10^* J) 


6 param. weight. 


1.6 


9.8 


1.91 


9.6680 


0.090 


-0.092 


890 


9.7 


6 param. unweight. 


1.6 


8.1 


2.07 


9.6677 


0.106 


-0.064 


860 


8.8 


Isothermal weight. 


1.5 


7.4 


2.38 


9.6679 


0.115 




860 


8.3 


Isothermal unweight. 


1.5 


11.0 


2.38 


9.6691 


0.078 




860 


8.4 


SN 1987A 


(10- kg) 


(10^^ m) 


(10-^ K) 


(MJD) 


(10-** s-^) 




(km s"-') 


(10*'' J) 


5-band UBVRI 


2.6 


2.1 


4.72 


47018.7 


7.2 


0.066 


1500 


4.5 


4-band BVRI 


2.5 


2.2 


4.69 


47025.3 


7.0 


0.071 


1500 


4.3 



Table 1: Derived and auxiliary parameters from the analytic fits to flickering of SS Cyg and the lightcurve 
of SN 1987A. 



Comparison of Model and Data Spectra 




4000 5000 6000 



Wavelength (A) 

Fig. 16. — Comparison of 3 spectra at UT 9.6722 (rise), 9.6830 (peak) and 9.6916 (fall) calculated from numerical 
integration across our fireballs with the observed spectra of SS Cyg. The '6 parameter weighted' values were 
used to generate the numerical spectra. 
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parameters rather than fixing the combination of 
to and H to give a launch time of JD 2446849.82 
when the associated neutrino event was detected 
by Kamiokande-II (Koshiba 1987). Since wc have 
not attempted to remove a "background" level the 
fit may be contaminated by the early stages of the 
radioactive tail. However, the flux in this region 
is dominated by the expanding/collapsing photo- 
sphere effect and any heating should show up in 
the fitted value of F differing from zero. 

We see in Figure 17 that we can achieve a qual- 
itative match between the analytic expression and 
the data. However, the U band is particularly 
poorly reproduced, representing the effect of blan- 
keting noted by Menzies et al. (1987). We refitted 
to the data using just the BVRI bands and show 
the results in Figure 18. The parameters for the 
fits are included in Table 5.1. We see that elimi- 
nating the U band data has relatively little effect 
and derive a value for the mass of the material in 
the expansion of M = 2.6 x 10^^ kg (« 13Mq). 
This compares to an estimate of 7Mq < M^nv < 
UMq by Saio, Nomoto & Kato (1988). The 
typical expansion velocity of 1500 km s~^ com- 
pares to a value derived from measurements of the 
blueshifts in P Cygni profiles of ^ 2500 km s~^ 
at maximum light (McCray 1993). Typical esti- 
mates of ii^kin ^ 2 X 10^^ J (Bethe & Pizzochero 
1990; Woosley 1988) are somewhat larger than 
our E^kin ~ 4 X lO''^ J. 

It is unsuprising that the lightcurves do not 
match in detail, since the Gaussian density pro- 
file we have assumed does not accurately repre- 
sent supernova ejecta. In particular, the reversal 
of the fast rise and slow decline derived analyt- 
ically can be understood from the way in which 
a photosphere marching inwards through a shell 
(in comoving coordinates) evolves. In the shell 
case, the material will first become optically thin 
at low impact parameters once the photosphere 
reaches the inner surface of the shell and later at 
higher impact parameters with their greater path- 
lengths through the shell material. While these 
can not be considered in any sense a good fits, we 
do retrieve values for the mass, size and temper- 
ature (see Figure 19) for this phase comparable 
with those given by Catchpole et al. (1987). We 
also see the same isothermal behaviour from the 
lightcurve analysis as those authors. The size of 
our photosphere is less than the radius they de- 



rive which is consistent with our including both 
optically thick and thin emission regions whereas 
they treat all the emission as optically thick. Our 
lightcurves also sometimes show an intrigiiing fea- 
ture (unfortunately most apparent in the poorly 
fitted U band) of a point of inflexion at around 
JD 2446960. This is reminiscent, in timing and 
in character, of the break in the lightcurves at- 
tributed to heating by radioactive decay. In our 
lightcurves, however, it appears to arise from the 
point at which the material becomes completely 
optically thin with no contribution from an opti- 
cally thick core region. All these factors, showing 
reasonable results for SN 1987A despite the inher- 
ent simplicity of our model, reinforce our confi- 
dence in the robustness of the parameters derived 
for SS Cyg which had much better agreement be- 
tween data and theory. 

6. Summary 

We have derived a formalism and analytic ex- 
pressions applicable to a variety of systems that re- 
produces the evolution of their optical lightcurves. 
The method has been tested for two widely dif- 
ferent cases allows us to derive reasonable values 
for the physical parameters involved in the ex- 
pansion. There is encouraging agreement with 
a method using a full integration of the opac- 
ity across the expanding region to generate both 
lightcurves and spectra. We have seen how the 
flickering of close binary systems can be modelled 
as smaller, hotter analogues of the well known 
rebrightening "bump" in supernovae lightcurves. 
We hope to test this approach further by compar- 
ison to data from LMXBs and other systems in 
the future. 

We thank the referee for the thoughtful and 
helpful suggestions in response to an earlier ver- 
sion of this paper. KJP has been supported, in 
part, by the U.S. National Science Foundation 
through grant AST-9987344 and, in part, through 
LSU's Center for Applied Information Technology 
and Learning. He also thanks Juhan Prank for 
stimulating discussions that improved this paper. 
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Fig. 8. — The value of ni at 5000 A for tempera- 
tures ranging in 2000 K steps from 8000 K (upper 
left curve) to 24 000 K. 
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Fig. 9. — Schematic representation of the alter- 
native integration schemes. The dot-dashed lines 
mark the two possible lines of sight considered for 
an observer (at x = 00). Along the y = line of 
sight the one-zone model would calculate (3^ us- 
ing «(= Kbf) evaluated at the photosphere. For 
a three zone scheme 0c would be calculated us- 
ing a combination of the relevant opacities evalu- 
ated at the indicated positions. The boundaries 
for each zone occur at the intersection of the dot- 
dashed line with the solid lines. Similarly the 
position where the opacities would be calculated 
for an optically thick/thin transition line of sight 
are indicated. While the circles of rj^ have been 
placed purely for illustrative convenience, the po- 
sition of the photosphere has been calculated as- 
suming To = 10. Other possible combination exist 
depending on the relative positions of the photo- 
sphere, peak density and ionization boundaries. 




Fig. 10. — Comparison of predicted lightcurves 
using the 4 combinations of 1- or 3-zone and 
y = or j/ni line of sight integration methods. 
The 4 methods produce virtually indistinguish- 
able results. The parameters of the model are 
M = S.OxlQi^ kg, ao = 7.5x10^ m, To = 16 000 K 
and r = 0. A lightcurve calculated from a full 
numerical integration of the opacity through the 
fireball is plotted as a dashed line. A distance of 
166 pc has been assumed. 
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Fig. 11. — The same comparison as Figure 10 but 
with To = 10 000 K.The 1- and 3- zone models are 
indistinguishable in either case. The y = models 
produce slightly lower peak fluxes. 
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Fig. 12.— Analytic fits (solid) to the SS Cyg data 
points and lightcurves generated from numerically 
calculated spectra using the best-fit parameters 
(dashed). F was allowed to be a fit parameter 
and the points were weighted according to the ob- 
servational errors. 
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Fig. 14. - Same comparison as Figmc 12, T = 
(isothermal) was fixed and the points were 
weighted according to the observational errors. 
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Fig. 13. — Same comparison as Figure 12, P was 
allowed to be a fit parameter and the points were 
weighted equally. 
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Figure 12, F = 



(isothermal) was fixed and the points were 
weighted equally. 



17 



M=.255E+32 kg ap=.213E-H4 r 
Time (JD 2446B0D+) 



Tg= 4719. K tQ=219.2 H = .723E-07 s ' r=0.0657 
"ime (JD 24468004) Time (JD 244680D+) 



u 

^HWiifiiiiii^ 1 i" ["iirViTnii 


B 




R 




\ /' 1 1 1 ..L^ 

* / \, 



Time (JD 2446800+) 



Time (JD 2446800+) 



Time (JD 2446800+) 



Fig. 17.— Analytic fits (solid) to the UBVRI 
SN 1987A data of Catchpole et al. (1987). T 
was allowed to be a fit parameter and all the dat- 
apoints with f > JD 2446875.0 were given equal 
weight. 
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Fig. 18. — Same comparison as Figure 17 fitting 
only BVRI band data. 
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Fig. 19. — The temperature and radius derived by 
Catchpole et al. (1987) (asterisks) and temper- 
ature, optically thick radius (in the V band) and 
current lengthscale from our UBVRI fit parame- 
ters. 
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